home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 2 / Apprentice-Release2.iso / Source Code / Libraries / VideoToolbox 94.11.17 / VideoToolboxSources / Luminance.c < prev    next >
Encoding:
Text File  |  1994-10-25  |  39.3 KB  |  1,065 lines  |  [TEXT/KAHL]
  1. /*
  2. Luminance.c
  3. See Luminance.h for prototypes.
  4. See Luminance.note for explanation.
  5. Also see:
  6. D.G. Pelli and L. Zhang (1991) Accurate control of contrast on microcomputer displays. 
  7. Vision Research, 31:1337-1360.
  8.  
  9. Copyright (c) 1989-1994 Denis G. Pelli. This software is free; you may use it in
  10. your research and give it away to others, with the following restrictions. Any
  11. copy you give away must include this paragraph, unmodified, and, if you have
  12. made any changes to the rest of the file, must include a note, added to HISTORY
  13. below, giving your name, the date, and a description of the changes. This
  14. software may not be sold, whether in source or compiled form, without my
  15. permission. I hope you will find this software useful, but I can't promise that
  16. it will work for you, and am not offering any support. That's why it's free. I
  17. would appreciate reports of bugs and improvements.
  18.  
  19. Denis G. Pelli, Ph.D.
  20. Institute for Sensory Research, Syracuse University, Syracuse, NY 13244-5290, U.S.A.
  21. denis_pelli@isr.syr.edu
  22.  
  23. HISTORY:
  24.  
  25. 4/24/89 dgp first working version.
  26. 4/28/89 dgp added device argument.
  27. 6/23/89 dgp added support for a video attenuator that has unequal gain in the three
  28.             channels. Added LoadClut() and LuminanceClut() routines.
  29. 7/30/89 dgp replaced call to SetEntries (which calls the Color Manager) by a call
  30.             to GDSetEntries (which calls the device driver).
  31. 8/10/89 dgp fixed bugs.
  32. 8/13/89 dgp rewrote LToV() to use nth order polynomial. Broke up LinearizeClut into
  33.             three routines--SetLuminanceRange, SetLuminance, and SetLuminances--the
  34.             last of which is equivalent to the old LinearizeClut.
  35. 8/16/89    dgp    After reading the Brooktree DAC manual I changed the documentation and
  36.             increased tolerance from 0.5 to 1 LSB, and now use the highest
  37.             luminance gain (worst case) over the requested range to transform it
  38.             into a luminance difference tolerance.
  39. 8/21/89    dgp    Made minor improvements in LToV() to deal with failure to converge due to 
  40.             getting stuck in a local minimum of the quartic. My solution is to now
  41.             require that Lmin and Lmax be filled in by the user in the calibration
  42.             structure. In practice this will guarantee that the minimum luminance
  43.             will not be below the local minimum found at the base of the rising curve.
  44.             Last week I also changed the LToV algorithm slightly, to do a bisection if
  45.             Newton's method would take us out of range.
  46. 9/7/89 dgp    Fixed bug in sorting routine, Sort().
  47. 9/10/89 dgp    Commented out the polynomial versions of VToL and LToV and replace them by
  48.             new versions that use a power law. The power law is a marginally better fit
  49.             than a 6th order polynomial (using only 4 instead of 7 parameters) and its
  50.             inverse can be computed much more quickly, about 0.3 ms instead of 2 ms.
  51. 9/11/89 dgp    Deciding to be indecisive, I introduced a conditional POWER_LAW_FIT
  52.             to control the choice of power law or polynomial fit for LToV and VToL.
  53. 9/16/89    dgp    Having finished writing the paper (Pelli & Zhang, 1991) I came back to
  54.             change the program to agree with what I wrote. I changed the equivalent
  55.             number tolerance to be the sum instead of the maximum of the 
  56.             variable-dac gains.
  57. 9/25/89 dgp Fixed minor bug in computation of tolerance. Now gives correct answer even
  58.             when the lowLuminance or highLuminance is out of range.
  59. 10/9/89 dgp    Made minor change so that when some of the gain[] are zero SetLuminance will
  60.             load zero into the unused colorSpec array elements. This is only a cosmetic
  61.             change.
  62. 10/28/89 dgp Made minor changes. I declared LToVPolynomial and LToVPower and
  63.             VToLPolynomial and VToLPower instead of having a conditional compilation. 
  64.             This makes both flavors of routine always available. I also changed 
  65.             LToVPolynomial to use LToVPower instead of LToVQuadratic for its initial 
  66.             guess. The quadratic guess was often poor and occasionally led to 
  67.             convergence problems. If the power law fit is not available then it reverts
  68.             to using the quadratic fit. If that's not available then it reverts to using
  69.             a middle-of-range guess.
  70. 10/28/89 dgp I eliminated the compile-time flag POWER_LAW_FIT and instead use whichever, 
  71.             polynomial or power, provides a better fit, as determined at run time. I 
  72.             biased the comparison of rms error to favor the power law fit since it has 
  73.             few parameters and can be inverted much more quickly. Note: if the 
  74.             power, polynomial, or quadratic fit parameters are not supplied then the 
  75.             appropriate LR.powerError, LR.polynomialError, or LR.quadraticError field 
  76.             should be set to infinity: INF or 1.0/0.0.
  77. 12/5/89    dgp    Added the trivial routine LToL() which enforces the bounds LP->LMin and 
  78.             LP->LMax.
  79. 4/2/90    dgp    Capitalized the word "To" in all function names, e.g. LToL().
  80. 4/22/90    dgp    Version 1.4. Made minor changes to the documentation. Renamed LToV() to
  81.             LToVFormulaic(), and introduced a new LToV() that uses a table to run faster.
  82.             LToV() takes an average of 170 µs, whereas LToVFormulaic() takes
  83.             1700 µs. I also now use the luminanceTable, if available, to speed up
  84.             VToL() slightly, from 310 µs to 180 µs. Each call to SetLuminance() now
  85.             takes 0.9 ms whereas it used to take 2.5 ms. The loss in accuracy is
  86.             negligible. The interpolation error can be reduced still further by
  87.             increasing LUMINANCES_IN_TABLE. Each doubling of the table size quarters the
  88.             maximum possible error of the interpolation.
  89. 4/23/90    dgp    Added a few lines of code to LToV() to check near the last index, in the
  90.             spirit of Numerical Recipes in C, hunt.c, before starting the bisection
  91.             search. I updated the times in the 4/22/90 note to reflect the latest timing
  92.             by TestLuminance.
  93. 7/27/90    dgp    Tightened up the error checking for illegal entries into the clut. I have
  94.             the impression that, contrary to specification, the Apple video driver
  95.             crashes and corrupts itself if given an out-of-range entry value in a
  96.             setEntry Control call.
  97. 9/18/90    dgp    Changed all instances of "v" to "V". The final version of the Pelli
  98.             & Zhang (1991) paper refers to the nominal voltage v; this file
  99.             now refers to the "equivalent number" V; they are related by V=255*v.
  100. 9/24/90    dgp    Updated the documentation to correspond more closely to the (hopefully)
  101.             final version of the Pelli & Zhang (1991) paper.
  102. 10/29/90 dgp Doubled speed of SetLuminance(). Now SetLuminance() takes 133 ms to
  103.             do a complete clut for a small new luminance range, and takes 188 ms for a
  104.             large new range. This required minor changes to SetLuminanceRange() & LToV().
  105.             Changed all function headers to Standard C prototype style.
  106. 10/30/90 dgp Replaced phrase "clut value" by "nominal voltage", for consistency with
  107.             Pelli & Zhang (1991).
  108.             Introduced new fixed fraction data type called "Milli", defined in Milli.h.
  109.             By judicious replacement of double by Milli, particularly in the luminance
  110.             table, I have increased the speed of SetLuminance by a factor of three.
  111.             A Mac II now takes only 42 to 50 ms to build a whole clut. This new feature
  112.             is enabled by setting FAST_LUMINANCE to 1 in Luminance.h. There is a slight
  113.             increase in the luminance tolerance.
  114. 10/31/90 dgp More fine tuning. Now takes 31 to 44 ms to build a whole clut. Introduced
  115.             conditional FAST_LUMINANCE in Luminance.h that causes the same code to
  116.             be compiled with either Milli or double variables. The _SetLuminance()
  117.             routine is now quite polished, and runs fast. All variables and routines
  118.             whose names begin with underscore _ are written here with type Milli, but
  119.             if FAST_LUMINANCE is false then this type is redefined (in this file only)
  120.             as double, and all the Milli macros are redefined to be appropriate for
  121.             a double. So bear in mind that things beginning with underscore are of
  122.             unknown type. The virtue of being able to run the same code with Milli or
  123.             double computations is that if you suspect an error in the (homemade)
  124.             Milli arithmetic, you can try out double arithmetic simply by setting
  125.             FAST_LUMINANCE to 0 and recompiling.
  126.             In the interest of speed I have streamlined the tolerance
  127.             computation. The new tolerances are probably as good as the old ones.
  128. 11/2/90 dgp Added _SetLuminances() subroutine that substitutes linear interpolation
  129.             for most of the calls to _LToV(). This has greatly speeded up
  130.             SetLuminancesAndRange(), with minor loss of accuracy. A complete clut now
  131.             takes 9 to 31 ms. Even a Mac II can now compute low-contrast cluts
  132.             on the fly, between frames. (The user may wish to optimize the speed-
  133.             accuracy trade-off controlled by the LINEAR_V_DOMAIN parameter, which
  134.             determines the maximum width of the interpolation interval.)
  135. 11/6/90 dgp Replaced Milli by FIXED, which can be compiled as either double or Fixed.
  136.             The Fixed math is sufficiently precise as to give the same DAC values as
  137.             double math.
  138.             Now figure out optimal bit shift LT->L.LShift to preserve as much resolution
  139.             as possible when interpolating in _LToV(). 
  140.             A complete clut now takes 8 to 30 ms. For threshold stimuli the average
  141.             time will usually be less than 15 ms.
  142. 11/8/90 dgp    Eliminated gamma slope table since the speed-up it offered was too small
  143.             to measure. Tidied up the computation of LShift.
  144. 8/24/91    dgp    Made compatible with THINK C 5.0.
  145. 12/17/91 dgp Added new routines: IncrementLuminance() and GetV().
  146.             The former is used to obtain the lowest possible contrast. On a log contrast
  147.             scale, minus infinity (zero contrast) is a poor approximation to even
  148.             a small log contrast.
  149. 3/11/92    dgp    Minor editing of comments.
  150. 12/2/92 dgp Fixed minor bug in SetLuminanceRange() reported by Wei Xie and Ken Alexander
  151.             that could cause the THINK C Debugger to report an error when an 
  152.             uninitialized double V[] was converted to an int. The value was not used 
  153.             for anything, so this was innocuous.
  154. 12/17/92 dgp 1.94 Began enhancement to allow any dacSize, not just 8-bit dacs. I have
  155.             confirmed that everything still works fine with 8-bit dacs, but I don't
  156.             have a 9-bit dac video card to test it any further. •CAUTION: this probably
  157.             does not yet work correctly for 9-bit dacs. E.g. it seems likely that the 
  158.             fixed point math will break if VMax is increased to 511; a lot of careful 
  159.             error analysis went into the original coding, at which time I thought that 
  160.             it was safe to assume an 8-bit dac (i.e. VMax==255). •Removed obsolete 
  161.             support for THINK C 4. 
  162. 12/21/92 dgp No longer load unused dac bits.
  163. 1/4/93    dgp    LoadLuminances now checks the gdType.
  164. 5/10/93    dgp    If driver is in gray mode (i.e. not color) then LoadLuminances maps rgb to 
  165.             gray by simply copying the green component to the red and blue components.
  166. 5/12/93    dgp removed debugging code that forced isGray to always be true.
  167. 3/24/94    dgp    added LToE, LToEOrdered, and EToL.
  168. 6/29/94 dgp renamed "LtoEOrdered" to "LToEOrdered", fixing erroneous capitalization.
  169. 9/28/94    dgp    If necessary, IncrementLuminance now bumps up the upper luminance bound to 
  170.             include the incremented luminance.
  171. 10/25/94 dgp LoadLuminance and GetV now ignore a device with no driver. (E.g. a device
  172.             created by NewGWorld.)
  173. */
  174. #include "VideoToolbox.h"        /* prototype for GDSetEntries() */
  175. #include "Luminance.h"
  176. #include <assert.h>
  177. #include <math.h>
  178.  
  179. static void Sort(int n,double arr[],int krr[]);    /* for internal use only */
  180.  
  181. #undef LongToFix
  182. #undef FixToLong
  183. #undef DoubleToFix
  184. #undef FixToDouble
  185. #if FAST_LUMINANCE
  186.     #define FIXED Fixed
  187.     #if !MACINTOSH
  188.         typedef long Fixed;
  189.         #define FixMul(x,y) DoubleToFix(FixToDouble(x)*FixToDouble(y))
  190.         #define FixDiv(x,y) DoubleToFix(FixToDouble(x)/FixToDouble(y))
  191.     #endif
  192.     #define LongToFix(x) ((long)(x)<<16)
  193.     #define FixToLong(x) ((x)>>16)
  194.     #define DoubleToFix(x) ((Fixed)((x)*65536.+0.5))
  195.     #define FixToDouble(x) ((double)(x)*(1./65536.))
  196.     #define D(x) FixToDouble(x)                            /* handy for debugging */
  197. #else
  198.     #define FIXED double
  199.     #define FixMul(x,y) ((double)(x)*(y))
  200.     #define FixDiv(x,y) ((double)(x)/(y))
  201.     #define LongToFix(x) ((double)(x))
  202.     #define FixToLong(x) ((long)(x))
  203.     #define DoubleToFix(x) ((double)(x))
  204.     #define FixToDouble(x) ((double)(x))
  205. #endif
  206.  
  207. #undef MAX
  208. #undef MIN
  209. #define MAX(a,b) ((a) > (b) ? (a) : (b))
  210. #define MIN(a,b) ((a) < (b) ? (a) : (b))
  211.  
  212. double EToL(LuminanceRecord *LP,int entry)
  213. {
  214.     return GetLuminance(NULL,LP,entry);
  215. }
  216.  
  217. int LToE(LuminanceRecord *LP,double L,int firstEntry,int lastEntry)
  218. //Returns the index of the table entry in specified range with luminance closest to L.
  219. {
  220.     int i,entry;
  221.     double dL,latestDL;
  222.     
  223.     latestDL=INF;
  224.     entry=0;
  225.     for(i=firstEntry;i<=lastEntry;i++){
  226.         dL=fabs(L-GetLuminance(NULL,LP,i));
  227.         if(dL<latestDL){
  228.             entry=i;
  229.             latestDL=dL;
  230.         }
  231.     }
  232.     return entry;
  233. }
  234.  
  235. int LToEOrdered(LuminanceRecord *LP,double L,int firstEntry,int lastEntry)
  236. /*
  237. Returns the index of the table entry in specified range with luminance closest to L.
  238. Assumes ordered table from firstEntry to lastEntry, either increasing or decreasing.
  239. Returned entry is guaranteed to be in the range firstEntry...lastEntry.
  240. Based on Numerical Recipes "locate.c", page 117
  241. */
  242. {
  243.     register int lower,middle,upper;
  244.     register double LLower,LMiddle,LUpper;
  245.     int ascend;
  246.     
  247.     lower=firstEntry;
  248.     upper=lastEntry;
  249.     LLower=GetLuminance(NULL,LP,lower);
  250.     LUpper=GetLuminance(NULL,LP,upper);
  251.     ascend=(LUpper>LLower);
  252.     while(upper-lower>1){
  253.         middle=(upper+lower)/2;
  254.         LMiddle=GetLuminance(NULL,LP,middle);
  255.         if(L>LMiddle == ascend){
  256.             lower=middle;
  257.             LLower=LMiddle;
  258.         }
  259.         else{
  260.             upper=middle;
  261.             LUpper=LMiddle;
  262.         }
  263.     }
  264.     if(L-LLower<LUpper-L == ascend)return lower;
  265.     else return upper;
  266. }
  267. double SetLuminance(GDHandle device,LuminanceRecord *LP
  268.     ,int theEntry,double luminance
  269.     ,double lowLuminance,double highLuminance)
  270. /*
  271. Set one entry in the ColorSpec table (and the clut if device is not NULL) to
  272. a specified luminance. It's ok for lowLuminance to be greater than highLuminance; 
  273. they still designate a range.
  274. */
  275. {
  276.     FIXED _luminance;
  277.     
  278.     if(theEntry<0 || theEntry>=COLORS 
  279.         || device!=NULL && theEntry>=GDCLUTSIZE(device))
  280.             PrintfExit("\007SetLuminance: illegal entry %d\n",theEntry);
  281.     
  282.     SetLuminanceRange(LP,lowLuminance,highLuminance);
  283.     _luminance=DoubleToFix(luminance);
  284.     _SetLuminance(LP,theEntry,_luminance);
  285.     if(device != NULL)LoadLuminances(device,LP,theEntry,theEntry);
  286.     return FixToDouble(_Tolerance(LP,_luminance));
  287. }
  288.  
  289. double SetLuminances(GDHandle device,LuminanceRecord *LP
  290.     ,int firstEntry,int lastEntry
  291.     ,double firstLuminance,double lastLuminance)
  292. /*
  293. Set a series of entries in the ColorSpec table (and the clut if device is not NULL)
  294. to a linear sequence of luminances. Assume this is the entire luminance range of 
  295. interest.
  296. */
  297. {
  298.     return SetLuminancesAndRange(device,LP,firstEntry,lastEntry
  299.         ,firstLuminance,lastLuminance,firstLuminance,lastLuminance);
  300. }
  301.  
  302. double SetLuminancesAndRange(GDHandle device,LuminanceRecord *LP
  303.     ,int firstEntry,int lastEntry
  304.     ,double firstLuminance,double lastLuminance
  305.     ,double lowLuminance,double highLuminance)
  306. /*
  307. Set a series of entries in the ColorSpec table (and the clut if device is not
  308. NULL) to a linear sequence of luminances. Uses last two arguments to set the
  309. luminance range.
  310.  
  311. I introduced a stack-space check before calling _SetLuminances(), since it calls
  312. itself recursively and may use 1000 bytes all together. However, printing an
  313. error message requires at least 4500 bytes. So I quit if the stack gets below
  314. 6000, so that the user will be presented with an understandable error message
  315. rather than a mysterious quit or crash.
  316. */
  317. {
  318.     FIXED _tolerance,_t,_dL64;
  319.     FIXED _firstL,_lastL,_firstV,_lastV;
  320.     
  321.     if(firstEntry<0 || firstEntry>lastEntry || lastEntry>=COLORS 
  322.         || device!=NULL && lastEntry>=GDCLUTSIZE(device) )
  323.         PrintfExit("\007SetLuminancesAndRange: illegal entries %d %d\n",firstEntry,lastEntry);
  324.     if(StackSpace()<6000)PrintfExit("SetLuminancesAndRange: not enough stack space.\007\n");
  325.     SetLuminanceRange(LP,lowLuminance,highLuminance);
  326.     _firstL=DoubleToFix(firstLuminance);
  327.     _lastL=DoubleToFix(lastLuminance);
  328.     _firstV=_LToV(LP,_firstL + LP->_LOffset);
  329.     _lastV=_LToV(LP,_lastL + LP->_LOffset);
  330.     if(lastEntry != firstEntry){
  331.         #if FAST_LUMINANCE
  332.         _dL64=DoubleToFix(64.0*(lastLuminance-firstLuminance)/(lastEntry-firstEntry));
  333.         #else
  334.         _dL64=DoubleToFix((lastLuminance-firstLuminance)/(lastEntry-firstEntry));
  335.         #endif
  336.     }else _dL64=0;
  337.     _SetLuminances(LP,firstEntry,lastEntry,_firstL,_dL64,_firstV,_lastV);
  338.     if(device != NULL)LoadLuminances(device,LP,firstEntry,lastEntry);    
  339.     
  340.     /* quick and dirty estimate of tolerance */
  341.     _tolerance=LP->L._Lu[0] - _firstL;
  342.     _t=_lastL - LP->L._Lu[LUMINANCES_IN_TABLE-1];
  343.     if(_tolerance<_t)_tolerance=_t;
  344.     if(_tolerance<0)_tolerance=0;
  345.     _tolerance+=LP->_tolerance;
  346.     return FixToDouble(_tolerance);    /* estimate of max error in ΔL */
  347. }
  348.  
  349. void _SetLuminances(LuminanceRecord *LPtr,int first,int last
  350.     ,FIXED _firstL,FIXED _dL64,FIXED _firstV,FIXED _lastV)
  351. /*
  352. This routine eliminates most of the the slow _LToV() table lookups and
  353. interpolations by assuming that the gamma function is smooth enough that we may
  354. linearly interpolate over any interval of V no wider than LINEAR_V_DOMAIN. Since
  355. the gamma function is roughly parabolic, the consequent error in L will be
  356. proportional to the square of LINEAR_V_DOMAIN, and independent of where this
  357. interval is along V. Since we're inscribing straight line segments in a
  358. positively accelerating gamma function, we will overestimate luminance, and
  359. consequently V will be too low. The error will be greatest at the middle of each
  360. linearly interpolated V interval.
  361.  
  362. If the requested V interval is larger than LINEAR_V_DOMAIN, then it is split
  363. into two intervals by making a pair of recursive calls. LINEAR_V_DOMAIN is
  364. defined in Luminance.h. I suggest a value of 4, but it may be set larger to
  365. attain slightly higher speed at lower accuracy.
  366.  
  367. CAUTION: This routine is a stack hog. Each call uses up about 64 bytes on the
  368. stack, and it calls itself recursively, up to log2(256/LINEAR_V_DOMAIN) times.
  369. Do NOT change any of the declarations of the "new" variables used in the first
  370. if{} to "static", as that would cause the recursion to fail.
  371. */
  372. {
  373.     static int doLast=1;
  374.     int saveDoLast;
  375.     int newOne;
  376.     FIXED _newL,_newV;
  377.     register LuminanceRecord *LP=LPtr;
  378.     register int i,Vi;
  379.     register FIXED _VToGo,_g;
  380.     static RGBColor *RGBPtr;    /* static in order to minimize stack usage */
  381.     static int theEntry;        /* static in order to minimize stack usage */
  382.     static FIXED _V,_dV;        /* static in order to minimize stack usage */
  383.     register short leftShift;
  384.     
  385.     if(last-first>1 ){
  386.         _dV=_lastV-_firstV;
  387.         if(_dV>LongToFix(LINEAR_V_DOMAIN) || _dV<LongToFix(-LINEAR_V_DOMAIN)){
  388.             newOne=(first+last)>>1;
  389.             #if FAST_LUMINANCE
  390.                 _newL=_firstL+((newOne-first)*_dL64>>6);
  391.             #else
  392.                 _newL=_firstL+(newOne-first)*_dL64;
  393.             #endif
  394.             _newV=_LToV(LP,_newL + LP->_LOffset);
  395.             saveDoLast=doLast;
  396.             doLast=0;
  397.             _SetLuminances(LP,first,newOne,_firstL,_dL64,_firstV,_newV);
  398.             doLast=saveDoLast;
  399.             _SetLuminances(LP,newOne,last,_newL,_dL64,_newV,_lastV);
  400.             return;
  401.         }
  402.     }
  403.     if(last!=first)_dV=(_lastV-_firstV)/(last-first);
  404.     else _dV=0;
  405.     if(!doLast)last--;
  406.     leftShift=LP->leftShift;
  407.     for(theEntry=first,_V=_firstV;theEntry<=last;theEntry++,_V+=_dV){
  408.         /****** This section of code is copied from _SetLuminance() ****/
  409.         RGBPtr = &LP->table[theEntry].rgb;
  410.         *RGBPtr=LP->rgb;                        /* load for fixed DACs */
  411.         _VToGo=_V - LP->_VFixed + LP->_VHalfStep;
  412.         for(i=LP->fixed;i<LP->dacs;i++) {
  413.             _g=LP->_gain[i];
  414.             Vi=_VToGo/_g;                        /* truncate to integer */
  415.             if(Vi>LP->VMax)Vi=LP->VMax;
  416.             if(Vi<LP->VMin)Vi=LP->VMin;
  417.             _VToGo -= Vi*_g;
  418.             ((short *)RGBPtr)[LP->dac[i]]=Vi<<leftShift;
  419.         }
  420.         /***************************************************************/
  421.     }
  422.     return;
  423. }
  424.  
  425. void _SetLuminance(LuminanceRecord *LPtr,int theEntry,FIXED _luminance)
  426. /*
  427. Set one entry in the ColorSpec table to a specified luminance. This is the
  428. private subroutine that actually does all the work for SetLuminance(). This
  429. routine is designed to run as fast as possible.
  430. */
  431. {
  432.     register LuminanceRecord *LP=LPtr;
  433.     register int i,Vi;
  434.     register FIXED _VToGo,_g;
  435.     RGBColor *RGBPtr;
  436.     register short leftShift=LP->leftShift;
  437.     
  438.     RGBPtr = &LP->table[theEntry].rgb;
  439.     *RGBPtr=LP->rgb;                        /* load for fixed DACs */
  440.     _VToGo=_LToV(LP,_luminance + LP->_LOffset) - LP->_VFixed + LP->_VHalfStep;
  441.     for(i=LP->fixed;i<LP->dacs;i++) {
  442.         _g=LP->_gain[i];
  443.         Vi=_VToGo/_g;                        /* truncate to integer */
  444.         if(Vi>LP->VMax)Vi=LP->VMax;
  445.         if(Vi<LP->VMin)Vi=LP->VMin;
  446.         _VToGo -= Vi*_g;
  447.         ((short *)RGBPtr)[LP->dac[i]]=Vi<<leftShift;
  448.     }
  449. }
  450.  
  451. void IncrementLuminance(GDHandle device,register LuminanceRecord *LP,int theEntry)
  452. /*
  453. Make smallest possible increase of the luminance of one entry in the ColorSpec
  454. table. If necessary, bump up the upper luminance bound to include the incremented luminance.
  455. */
  456. {
  457.     double V,L,highL;
  458.     
  459.     V=GetV(device,LP,theEntry);
  460.     V+=LP->VHalfStep;
  461.     V+=LP->VHalfStep;
  462.     L=VToL(LP,V);
  463.     if(L>LP->highLuminance)highL=L;
  464.     else highL=LP->highLuminance;
  465.     SetLuminance(device,LP,theEntry,L,LP->lowLuminance,highL);
  466. }
  467.  
  468. FIXED _Tolerance(LuminanceRecord *LPtr,FIXED _luminance)
  469. {
  470.     register LuminanceRecord *LP=LPtr;
  471.     register FIXED _tolerance,_t;
  472.     
  473.     /* quick and dirty estimate of tolerance */
  474.     if(LP->L.exists==luminanceSet){
  475.         _tolerance=LP->L._Lu[0] - _luminance;
  476.         _t=_luminance - LP->L._Lu[LUMINANCES_IN_TABLE-1];
  477.         if(_tolerance<_t)_tolerance=_t;
  478.         if(_tolerance<0)_tolerance=0;
  479.         _tolerance+=LP->_tolerance;
  480.     }
  481.     else _tolerance=LP->_tolerance;
  482.     return _tolerance;
  483. }
  484.  
  485. double GetLuminance(GDHandle device,LuminanceRecord *LP,int theEntry)
  486. /*
  487. If device is not NULL then examines one entry in the actual clut, otherwise
  488. examines the ColorSpec table contained in *LP, and in either case returns the
  489. luminance that will be produced. Supplying an illegal entry value results in a
  490. returned value of -INF.
  491. */
  492. {
  493.     return VToL(LP,GetV(device,LP,theEntry));
  494. }
  495.  
  496. double GetV(GDHandle device,LuminanceRecord *LP,int theEntry)
  497. /*
  498. If device is not NULL (and has a driver) then examines one entry in the actual clut, otherwise
  499. examines the ColorSpec table contained in *LP, and in either case returns the V
  500. that will be produced. Supplying an illegal entry value results in a returned
  501. value of -INF.
  502. */
  503. {
  504.     RGBColor *myRGBPtr=NULL;
  505.     ColorSpec myCSpec;
  506.     double V;
  507.     int error;
  508.     
  509.     if(device != NULL && (*device)->gdRefNum!=0) {
  510.         if(theEntry<0 || theEntry>=GDCLUTSIZE(device)){
  511.             PrintfExit("GetLuminance: illegal entry %d\n",theEntry);
  512.         }
  513.         error=GDGetEntries(device,theEntry,0,&myCSpec);
  514.         if(error){
  515.             PrintfExit("GetLuminance: GDGetEntries error %d\007",error);
  516.         }
  517.         myRGBPtr=&myCSpec.rgb;
  518.     }else{
  519.         if(LP!=NULL && theEntry>=0 && theEntry<COLORS) 
  520.             myRGBPtr=&LP->table[theEntry].rgb;
  521.     }
  522.     if(myRGBPtr == NULL) return -INF;
  523.     V = LP->r*(myRGBPtr->red>>8);
  524.     V += LP->g*(myRGBPtr->green>>8);
  525.     V += LP->b*(myRGBPtr->blue>>8);
  526.     return V;
  527. }
  528.  
  529. void LoadLuminances(GDHandle device, LuminanceRecord *LP,
  530.     int firstEntry, int lastEntry)
  531. /*
  532. This just calls GDSetEntries() or GDDirectSetEntries() to load your ColorSpec
  533. table into the clut of your screen device. It is here simply to provide a
  534. cosmetic match to the call to SetLuminances(), for which loading the clut is
  535. optional. Note: if you prefer, instead of LP you may send just the address of a
  536. ColorSpec table, cast to (LuminanceRecord *), since a LuminanceRecord begins
  537. with a ColorSpec table.
  538. */
  539. {
  540.     int error;
  541.     short isGray,i;
  542.     RGBColor *RGBPtr;
  543.     ColorSpec table[256],*tablePtr;
  544.  
  545.     if(device==NULL || (*device)->gdRefNum==0)return;    // no device, or no device driver
  546.     isGray=!TestDeviceAttribute(device,gdDevType);
  547.     if(isGray){
  548.         for(i=firstEntry;i<=lastEntry;i++){
  549.             RGBPtr=&table[i].rgb;
  550.             RGBPtr->red=RGBPtr->green=RGBPtr->blue=LP->table[i].rgb.green;
  551.         }
  552.         tablePtr=&table[firstEntry];
  553.     }else tablePtr=&LP->table[firstEntry];
  554.  
  555.     switch((*device)->gdType){
  556.     case fixedType:
  557.         printf("LoadLuminances: this device has a fixed CLUT.\n");
  558.         return;
  559.     case clutType:
  560.         error=GDSetEntries(device,firstEntry,lastEntry-firstEntry,tablePtr);
  561.         break;
  562.     case directType:
  563.         error=GDDirectSetEntries(device,firstEntry,lastEntry-firstEntry,tablePtr);
  564.         break;
  565.     }
  566.     if(error) printf("\007LoadLuminances: error %d\n",error);
  567. }
  568.  
  569. double LToL(LuminanceRecord *LP,double L)
  570. {
  571.     if(L > LP->LMax) return LP->LMax;
  572.     if(L < LP->LMin) return LP->LMin;
  573.     return L;
  574. }
  575.  
  576. double VToL(LuminanceRecord *LP,double V)
  577. /*
  578. Return the luminance that would result from a nominal voltage.
  579. */
  580. {
  581.     return FixToDouble(_VToL(LP,DoubleToFix(V)));
  582. }
  583.  
  584. FIXED _VToL(LuminanceRecord *LP,FIXED _V)
  585. /*
  586. Return the luminance that would result from a nominal voltage.
  587. */
  588. {
  589.     register int i;
  590.     register luminanceTable *LT;
  591.     register FIXED _di;
  592.     double LF,VF;
  593.     
  594.     LT=&LP->L;
  595.     if(LT->exists==luminanceSet){
  596.         _di=FixDiv(_V-LT->_VMin,LT->_dV);
  597.         i=FixToLong(_di);
  598.         _di -= LongToFix(i);
  599.         if(i<0)return LT->_Lu[0];
  600.         if(i>=LUMINANCES_IN_TABLE-1)return LT->_Lu[LUMINANCES_IN_TABLE-1];
  601.         return LT->_Lu[i]+FixMul(_di,LT->_Lu[i+1]-LT->_Lu[i]);
  602.     }
  603.     VF=FixToDouble(_V);
  604.     if(LP->powerError < 2.0*LP->polynomialError) LF=VToLPower(LP,VF);
  605.     else LF=VToLPolynomial(LP,VF);
  606.     return DoubleToFix(LF);
  607. }
  608.  
  609. double LToV(LuminanceRecord *LP,double L)
  610. {
  611.     return FixToDouble(_LToV(LP,DoubleToFix(L)));
  612. }
  613.  
  614. FIXED _LToV(LuminanceRecord *LP,FIXED _Lu)
  615. /*
  616. New, faster version uses a tabulated version of the gamma function VToL(). This
  617. replaces the old function LToVFormulaic(), formerly called LToV(). A subtlety is
  618. that I tabulate the gamma function _Lu[]=VToL() rather than its inverse, because
  619. this yields a much lower upper bound on the error of the linearly interpolated
  620. L. This is true because the gamma function is roughly parabolic. Equally spaced
  621. samples of a parabolic function yield equal maximum error of linear
  622. interpolation in all intervals. The cost of tabulating _Lu[] rather than its
  623. inverse is that we must search the table in order to find the largest _Lu[] less
  624. than or equal to L. The search time is minimized by using a search procedure
  625. that begins the search at the place found by the last search, since successive
  626. calls to LToV() are likely to request similar luminances.
  627. */
  628. {
  629.     register luminanceTable *LT;
  630.     register int lo,hi,i;
  631.     register FIXED _dL,_dLTemp;
  632.     FIXED _V,_LLo;
  633.     
  634.     LT=&LP->L;
  635.     if(LT->exists != luminanceSet){
  636.         /* get domain of the monotonic section of the gamma function */
  637.         LT->_VMin=DoubleToFix(LToVFormulaic(LP,LP->LMin));
  638.         LT->_VMax=DoubleToFix(LToVFormulaic(LP,LP->LMax));
  639.         LT->_dV=(LT->_VMax-LT->_VMin)/(LUMINANCES_IN_TABLE-1);
  640.         _V=LT->_VMin;
  641.         for(i=0;i<LUMINANCES_IN_TABLE;i++,_V+=LT->_dV) LT->_Lu[i]=_VToL(LP,_V);
  642.         /* just in case, impose monotonicity, from center on out */
  643.         for(i=LUMINANCES_IN_TABLE/2;i>0;i--)
  644.             if(LT->_Lu[i-1] > LT->_Lu[i]) LT->_Lu[i-1]=LT->_Lu[i];
  645.         for(i=1+LUMINANCES_IN_TABLE/2;i<LUMINANCES_IN_TABLE;i++)
  646.             if(LT->_Lu[i-1] > LT->_Lu[i]) LT->_Lu[i]=LT->_Lu[i-1];
  647.         #if FAST_LUMINANCE    /* compute the slope dL/dV */
  648.             /* find shift that maximizes ΔL precision without overflow */
  649.             _dL=0;
  650.             for(i=0;i<LUMINANCES_IN_TABLE-1;i++){
  651.                 _dLTemp=LT->_Lu[i+1]-LT->_Lu[i];
  652.                 if(_dL<_dLTemp)_dL=_dLTemp;
  653.             }
  654.             LT->LShift=Log2L(LONG_MAX/_dL);
  655.         #endif                
  656.         LT->exists=luminanceSet;
  657.         LT->latestIndex=-2;    /* invalid, so hunt will start from scratch */
  658.     }
  659.     /* hunt for L in LT->_Lu[] */
  660.     /* first check at latestIndex, latestIndex±1, ... */
  661.     lo=-1;
  662.     hi=LUMINANCES_IN_TABLE;
  663.     i=LT->latestIndex;
  664.     if(i>lo && i<hi){
  665.         if(_Lu>LT->_Lu[i])
  666.             {lo=i;i++;}
  667.         else
  668.             {hi=i;i--;}
  669.     }
  670.     /* simple bisection search, see Numerical Recipes in C, pages 98-99. */
  671.     while(hi-lo>1){
  672.         i=(hi+lo)>>1;
  673.         if(_Lu>LT->_Lu[i])lo=i;
  674.         else hi=i;
  675.     }
  676.     LT->latestIndex=lo;
  677.     if(lo<0){
  678.         LT->latestIndex=0;            /* nearest legal index */
  679.         return LT->_VMin;
  680.     }
  681.     if(lo>=LUMINANCES_IN_TABLE-1)return LT->_VMax;
  682.     _LLo=LT->_Lu[lo];
  683.     _dL=_Lu-_LLo;
  684.     #if FAST_LUMINANCE
  685.         return LT->_VMin+lo*LT->_dV
  686.             +(_dL<<LT->LShift)/((LT->_Lu[lo+1]-_LLo<<LT->LShift)/LT->_dV);
  687.     #else
  688.         return LT->_VMin+LT->_dV*(lo+_dL/(LT->_Lu[lo+1]-_LLo));
  689.     #endif
  690. }
  691.  
  692. double LToVFormulaic(LuminanceRecord *LP,double L)
  693. /*
  694. Find a nominal voltage that would give luminance L. Failing that, it returns the
  695. closest possible value. The answer is computed by inverting the formula for the
  696. gamma function.
  697. */
  698. {
  699.     if(LP->powerError < 2.0*LP->polynomialError) return LToVPower(LP,L);
  700.     else return LToVPolynomial(LP,L);
  701. }
  702.  
  703. double VToLPolynomial(LuminanceRecord *LP,double V)
  704. /*
  705. Return the luminance that would result from a nominal voltage. Uses a polynomial
  706. equation L = p[0] + p[1]*V^1 +... of any degree. dgp 8/11/89
  707. */
  708. {
  709.     register double L,VV;
  710.     register int i,m;
  711.     
  712.     m=LP->coefficients;
  713.     L=0.0;
  714.     VV=1.0;
  715.     for(i=0;i<m;i++){
  716.         L+=LP->p[i]*VV;
  717.         VV*=V;
  718.     }
  719.     return L;
  720. }
  721.  
  722. double VToLPower(LuminanceRecord *LP,double V)
  723. /*
  724. Return the luminance that would result from a nominal voltage.
  725. Uses a rectified power law:
  726.  
  727.     L(V)=power[0]+Rectify(power[1]+power[2]*V)^power[3]
  728.     
  729.     Rectify(x)=x if x≥0
  730.     Rectify(x)=0 if x<0
  731.  
  732. dgp 10/28/89
  733. */
  734. {
  735.     double L;
  736.     
  737.     L=LP->power[1]+LP->power[2]*V;
  738.     if(L>0.0) L=LP->power[0]+pow(L,LP->power[3]);
  739.     else L=LP->power[0];
  740.     return L;
  741. }
  742.  
  743. double LToVPolynomial(LuminanceRecord *LP,double L)
  744. /*
  745. Find a nominal voltage that would give luminance L. Failing that, it returns the
  746. closest possible value. Solves a polynomial equation L = p[0] + p[1]*V^1 +... of
  747. any degree. First we find a quick approximate answer by solving the power law
  748. fit, LToVPower(). Then we use Newton's method to home in quickly on the solution
  749. to the polynomial fit. Newton's method is taken from Numerical Recipes in C. For
  750. a polynomial of degree 8 (which is what I recommend) this routine now does about
  751. 3? iterations and takes about 1.1? ms. The error in V is less than TOLERANCE,
  752. i.e. relative to full scale of 256 we get an accuracy of 14.6 bits. You can
  753. increase the TOLERANCE to reduce the number of iterations to make it slightly
  754. faster. A TOLERANCE of 1.0 (for 8-bit accuracy) takes 0.8? ms, not much of a
  755. savings in time, and a large loss in accuracy. Reducing TOLERANCE to 1e-8 yields
  756. an accuracy of nearly 35 bits and takes 1.5? ms, only slightly longer.
  757.  
  758. This routine will be called essentially once per clut entry, so computing a
  759. whole new clut with entries 0..255 will take 256*1.1? ms = 282? ms. Some
  760. experiments use only part of the clut, and will therefore take less time. In
  761. some of the experiments that do use the whole clut it may be desirable to
  762. compute the clut table in advance.
  763. dgp 8/11/89
  764. */
  765. {
  766.     #define TOLERANCE 0.01                    /* desired accuracy of V, see above */
  767.     #define JMAX 20
  768.     double a[MAX_COEFFICIENTS];
  769.     register int i;
  770.     int m,j;
  771.     register double f,df,u,VV,V,dV;
  772.     
  773.     if(L>LP->LMax) return LP->VMax;
  774.     if(L<LP->LMin) return LP->VMin;
  775.  
  776.     m=LP->coefficients;
  777.     if(m>MAX_COEFFICIENTS || m<1)
  778.         PrintfExit("LToVPolynomial: %d coefficients is too many or too few\007\n",m);
  779.     for(i=0;i<m;i++) a[i]=LP->p[i];
  780.     a[0]-=L;
  781.     if(LP->powerError < 0.2*LP->LMax) V=LToVPower(LP,L);    /* very good initial guess */
  782.     else {
  783.         if(LP->quadraticError < 0.2*LP->LMax) V=LToVQuadratic(LP,L);/* fair initial guess */
  784.         else PrintfExit("LToVPolynomial: neither power nor quadratic fit is available\007\n");
  785.     }
  786.     for(j=0;j<JMAX;j++){
  787.         /* evaluate function, i.e. the polynomial, and its derivative */
  788.         f=a[0];
  789.         df=0.0;
  790.         VV=1.0;
  791.         for(i=1;i<m;i++){
  792.             u=a[i]*VV;
  793.             df+=i*u;
  794.             f+=V*u;
  795.             VV*=V;
  796.         }
  797.         dV=f/df;    /* apply Newton's method */
  798.         if(V-dV<LP->VMin) dV=(V-LP->VMin)/2.0;    /* bisect */
  799.         if(V-dV>LP->VMax) dV=(V-LP->VMax)/2.0;    /* bisect */
  800.         V -= dV;
  801.         if(fabs(dV) < TOLERANCE) return V;
  802.     }
  803.     printf("LToVPolynomial(%7.2f): Warning, too many iterations. VToL(%5.1f)=%7.2f\n",L,V,VToL(LP,V));
  804.     return V;
  805.     #undef TOLERANCE
  806.     #undef JMAX
  807. }
  808.  
  809. double LToVPower(LuminanceRecord *LP,double L)
  810. /*
  811. Find a nominal voltage that would give luminance L. Failing that, it returns the
  812. closest possible value. Solves the rectified power law:
  813.     L(V)=power[0]+Rectify(power[1]+power[2]*V)^power[3]
  814.     Rectify(x)=x if x≥0
  815.     Rectify(x)=0 if x<0
  816. LToVPower takes about 300 microseconds on a Mac II with 8881 arithmetic. This
  817. routine will be called essentially once per clut entry, so computing a whole new
  818. clut with entries 0..255 will take 256*0.3 ms = 77 ms. Some experiments use only
  819. part of the clut, and will therefore take less time. In experiments that do use
  820. the whole clut it may be desirable to compute the clut table in advance.
  821. dgp 10/28/89
  822. */
  823. {
  824.     double V;
  825.     
  826.     if(L<LP->power[0]) L=LP->power[0];
  827.     V=(pow(L-LP->power[0],1.0/LP->power[3])-LP->power[1])/LP->power[2];
  828.     if(V>LP->VMax) V=LP->VMax;
  829.     return V;
  830. }
  831.  
  832. double LToVQuadratic(LuminanceRecord *LP,double L)
  833. /*
  834. This is a quick, approximate version of LToV() that is used by LToVPolynomial
  835. for an initial guess if no power law fit is available. Find nominal voltage that
  836. would give luminance L. Since a quadratic fit to the luminance calibration is
  837. only a fair approximation, this routine serves only to find a quick approximate
  838. answer to the root of a higher order polynomial fit. Solves the quadratic
  839. equation by method recommended in Numerical Recipes in C. My choice of root
  840. assumes that L is an increasing function of V over the operating range of the
  841. display. If concave (i.e. q[2]<0) use the larger root. If convex (i.e. q[2]>0)
  842. use the smaller root. In plain English, if the parabola is right side up then we
  843. take the right hand branch. If the parabola is upside down then we take the left
  844. hand branch. dgp 8/12/89
  845. */
  846. {
  847.     register double a,b,c,q;
  848.     double V;
  849.  
  850.     if(L>LP->LMax) return LP->VMax;
  851.     if(L<LP->LMin) return LP->VMin;
  852.     c=LP->q[0]-L;
  853.     b=LP->q[1];
  854.     a=LP->q[2];
  855.     if(a == 0.0) return -c/b;
  856.     if(b<0.0)
  857.         q=-0.5*(b-sqrt(b*b-4.0*a*c));
  858.     else
  859.         q=-0.5*(b+sqrt(b*b-4.0*a*c));
  860.     if(a*b<0.0)
  861.         V=q/a;
  862.     else
  863.         V=c/q;
  864.     return V;
  865. }
  866.  
  867. double SetLuminanceRange(LuminanceRecord *LPtr,double lowLuminance,double highLuminance)
  868. /*
  869. The user will never call this routine directly, as it is called automatically by
  870. SetLuminance and SetLuminances. This routine does the hard work of figuring out
  871. which DACs to fix to optimally represent a given luminance range. The goal is to
  872. fix the coarser DACs and vary the finer DACs to represent the image. The CHANGES
  873. in luminance (i.e. the contrast) will have the accuracy of the coarsest of the
  874. variable DACs. To avoid the overhead of figuring all this out every single time
  875. you call SetLuminance(), the results are stored temporarily in your
  876. LuminanceRecord, and are only recomputed if you request a different luminance
  877. range.
  878.  
  879. This function checks and returns immediately if you call it with the same
  880. luminance range that it already has stored in your LuminanceRecord.
  881. */
  882. {
  883.     register LuminanceRecord *LP=LPtr;
  884.     register int i;
  885.     double dL,dL2,LOffset;
  886.     double VSum,VFixed;
  887.     int dacs,fixed,dac[DACS];
  888.     double gain[DACS],V[DACS];
  889.     double VGoal;
  890.     double lowV,highV;
  891.     double tolerance;
  892.     double VHalfStep;
  893.     double dV;
  894.     unsigned short *RGBArray;
  895.     
  896.     /* -1. If necessary, swap low & high */
  897.     if(lowLuminance > highLuminance){
  898.         dL=lowLuminance;
  899.         lowLuminance=highLuminance;
  900.         highLuminance=dL;
  901.     }
  902.  
  903.     /* 0. Return right away if our work has already been done */
  904.     if(LP->rangeSet == luminanceSet &&
  905.         LP->lowLuminance==lowLuminance &&
  906.         LP->highLuminance==highLuminance)
  907.             return LP->tolerance;
  908.  
  909.     /* 1. sort the DACs by gain, discard any with zero or negative gain.
  910.     Gain is defined as the change in V when a single DAC is increased from 0 to 1.
  911.     Sort the gains so that gain[0]≥gain[1]≥gain[2]...
  912.     Note this order is opposite to that used in Pelli & Zhang (1991). The answers
  913.     are the same, but the notation is different.
  914.     */
  915.     if(LP->rangeSet != luminanceSet){
  916.         /* This only needs to be done once, even if the range changes. */
  917.         gain[0]=LP->r;
  918.         gain[1]=LP->g;
  919.         gain[2]=LP->b;
  920.         #if DACS>3
  921.             #error "Current implementation doesn't support more than 3 DACs."
  922.         #endif
  923.         for(i=0;i<DACS;i++)dac[i]=i;
  924.         Sort(DACS,gain,dac);            /* sort so that gain[0]≥gain[1]≥gain[2] */
  925.         dacs=0;
  926.         for(i=0;i<DACS;i++) if(gain[i]>0.0)dacs++;
  927.         VHalfStep=gain[dacs-1]/2.0;        /* half the smallest step */
  928.         for(i=0;i<COLORS;i++)LP->table[i].value=i;
  929.         LP->dacSize=Log2L(LP->VMax*2-1);    // in bits, rounded up
  930.         if(LP->dacSize!=8)printf(
  931.             "Caution: Your video card seems to have %ld-bit video dacs, and Luminance.c\n"
  932.             "has only been tested with 8-bit dacs.\n",LP->dacSize);
  933.         /* Apple's convention is to provide a 16-bit value to the dac. Luminance.c
  934.         saves time and space by computing values with only as much precision as
  935.         needed for the the video card's dac (i.e. 0 to VMax, or dacSize bits). This
  936.         value must be shifted left by leftShift bits to yield a standard 16-bit value.
  937.         Note that Apple normally recommends scaling the value to fill the entire
  938.         16-bit range, e.g. multiplying an 8-bit value by 0x101, but that doesn't
  939.         seem appropriate here, where we know the dac size, and are more concerned with
  940.         step size that with range.
  941.         */
  942.         LP->leftShift=16-LP->dacSize;
  943.     }
  944.     else {
  945.         for(i=0;i<DACS;i++){
  946.             gain[i]=LP->gain[i];
  947.             dac[i]=LP->dac[i];
  948.         }
  949.         dacs=LP->dacs;
  950.         VHalfStep=LP->VHalfStep;
  951.     }
  952.     for(i=0;i<DACS;i++)V[i]=0.0;    /* initialize fixed dac voltages */
  953.     
  954.     /* 2. Transform lowLuminance and highLuminance to lowV and highV. */
  955.     lowV=LToV(LP,lowLuminance);        /* takes 0.15 ms */
  956.     highV=LToV(LP,highLuminance);    /* takes 0.15 ms */
  957.  
  958.     /* 3. Designate the finest dacs as variable until we have enough to cover
  959.     the range, then designate the rest as fixed. (DACs 0..fixed-1 will be fixed.)
  960.     */
  961.     VGoal=fabs(highV-lowV);
  962.     VSum=0.0;
  963.     fixed=0;
  964.     for(i=dacs-1;i>=0;i--) {
  965.         if(VSum >= VGoal)fixed++;
  966.         VSum += gain[i]*(LP->VMax - LP->VMin);
  967.     }
  968.     /* sum gains of all variable dacs. This is called gVary by Pelli & Zhang (1991) */
  969.     tolerance=0.0;
  970.     for(i=dacs-1;i>=fixed;i--)tolerance+=gain[i];
  971.     /* scale by highest luminance gain in the requested range */
  972.     dL=fabs(VToL(LP,lowV+1.0)-VToL(LP,lowV));
  973.     dL2=fabs(VToL(LP,highV)-VToL(LP,highV-1.0));
  974.     tolerance *= MAX(dL,dL2);
  975.  
  976.     /* 4. Temporarily set the variable DACs to their midpoints. Now set the fixed DACs
  977.     to most accurately represent (lowV+highV)/2.
  978.     */
  979.     VGoal=(lowV+highV)/2.0;
  980.     VSum=0.0;
  981.     VFixed=(LP->VMax + LP->VMin)/2.0;
  982.     for(i=fixed;i<dacs;i++) VSum+=gain[i]*VFixed;    /* The mid-voltage of the variable DACs */
  983.     for(i=0;i<fixed;i++) {
  984.         V[i]=floor((VHalfStep+VGoal-VSum)/gain[i]);    /* the unattenuated voltage of DAC i */
  985.         V[i]=MAX(LP->VMin,MIN(LP->VMax,V[i]));
  986.         VSum += V[i]*gain[i];
  987.     }
  988.     VSum=0.0;
  989.     for(i=0;i<fixed;i++) VSum+=gain[i]*V[i];
  990.     VFixed=VSum;
  991.     
  992.     /* 5. The limited precision of the fixed DACs may result in a 
  993.     small offset between the requested and now available range. The offset will be
  994.     at most half a step of the finest of the fixed DACs. It seems reasonable
  995.     to offset the requested luminance range by up to that small amount to better fit
  996.     the now available range. */
  997.     LOffset=0.0;
  998.     if(fixed>0){
  999.         VSum=VFixed;
  1000.         for(i=fixed;i<dacs;i++) VSum += LP->VMin*gain[i];
  1001.         dV= VSum - MIN(lowV,highV);
  1002.         dV= MIN(dV,gain[fixed-1]/2.0);
  1003.         if(dV>0.0) LOffset+=VToL(LP,VSum+dV)-VToL(LP,VSum);
  1004.         VSum=VFixed;
  1005.         for(i=fixed;i<dacs;i++) VSum += LP->VMax*gain[i];
  1006.         dV= MAX(lowV,highV) - VSum;
  1007.         dV= MIN(dV,gain[fixed-1]/2.0);
  1008.         if(dV>0.0) LOffset-=VToL(LP,VSum+dV)-VToL(LP,VSum);
  1009.     }
  1010.     /* Now save this information in the LuminanceRecord for future use */
  1011.     for(i=0;i<DACS;i++){
  1012.         LP->dac[i]=dac[i];
  1013.         LP->gain[i]=gain[i];
  1014.     }
  1015.     LP->rangeSet=luminanceSet;
  1016.     LP->lowLuminance=lowLuminance;
  1017.     LP->highLuminance=highLuminance;
  1018.     LP->dacs=dacs;
  1019.     LP->VHalfStep=VHalfStep;
  1020.     LP->fixed=fixed;
  1021.     LP->VFixed=VFixed;
  1022.     LP->tolerance=tolerance;
  1023.     LP->LOffset=LOffset;
  1024.  
  1025.     /* Copy into FIXED variables */
  1026.     for(i=0;i<DACS;i++)LP->_gain[i]=DoubleToFix(LP->gain[i]);
  1027.     LP->_VHalfStep=DoubleToFix(VHalfStep);
  1028.     LP->_VFixed=DoubleToFix(VFixed);
  1029.     LP->_tolerance=DoubleToFix(tolerance);
  1030.     LP->_LOffset=DoubleToFix(LOffset);
  1031.  
  1032.     /* cache the fixed DAC values, and zero the rest for good measure */
  1033.     RGBArray = (unsigned short *) &LP->rgb;    /* treat as an array */
  1034.     for(i=0;i<DACS;i++){
  1035.         RGBArray[LP->dac[i]] = (short)V[i]<<LP->leftShift;
  1036.     }
  1037.     return tolerance;
  1038. }
  1039.  
  1040. static void Sort(int n,double arr[],int krr[])
  1041. /*
  1042. Slightly modified sort routine piksr2() from Numerical Recipes in C. I changed
  1043. "float" to "double", and I changed second array to int. I reversed the order, so
  1044. largest element would be first. I changed it to use conventional c arrays,
  1045. starting at index 0.
  1046. */
  1047. {
  1048.     register int i,j;
  1049.     double a;
  1050.     int k;
  1051.  
  1052.     for(j=1;j<n;j++) {
  1053.         a=arr[j];
  1054.         k=krr[j];
  1055.         i=j-1;
  1056.         while (i >= 0 && arr[i] < a) {
  1057.             arr[i+1]=arr[i];
  1058.             krr[i+1]=krr[i];
  1059.             i--;
  1060.         }
  1061.         arr[i+1]=a;
  1062.         krr[i+1]=k;
  1063.     }
  1064. }
  1065.